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'—'"'ktr  r,r>  ^  A-,  ,  j  1  I  £  rsadJust  the  ••alman  «aln  (Ref-  11  •  However  this 

BSTRACTj  M  F —  I  /  m  Oz  ;  *  yT7  approach  was  not  very  fruitful  for  the  following 

I  '  _  * _ :  reasons:  Conventional  one-dimensional  (1-D)  re- 

■Vhen  an  image  is  estimated  from  noisy  data  \  /a — .  c’ursive  least  squares  identification  algorithms  do 
sing  a  linear  shift-invariant  (LSI)  filter,  the  j  If H  j  hot  readily  extend  to  2-D  recursive  identification 
ubjective  improvement  is  relatively  poor  at  low  ;'7  /algorithms.  Secondly,  in  order  that  the  recursive 
ignal-to-noise  ratios.  This  occurs  for  at  least  *  identification  algorithms  he  able  to  track  space 
.vo  reasons:  first,  the  statistics  of  the  image  varying  parameters,  the  effective  bandwidth  must  be 

o  are  markedly  space-variant  and  second,  the  eye  is  large  enough  that  the  variations  are  not  smoothed 

very  sensitive  to  blurring  of  edges.  However  adap-  over.  This  in  turn  makes  the  algorithm  more  sus- 

'‘tive  filtering  techniques  can  be  applied  to  improve  ceptible  to  noise  effects.  Thirdlv,  the  existence 

the  subjective  quality  of  noisy  images  even  at  low  of  process  noise  with  an  'unknown  v  iance  and  meas- 
signal-to-noise  ratios.  This  is  accomplished  in  urement  noise  as  well,  can  result  in  biased  estimates 

the  present  work  by  using  multiple  models  to  match  ir  General  least  squares  regression  is  performed 

the  space-variant  statistics  and  by  using  oriented  directly  on  the  measurements, 

edge  models  to  prevent  edge  blurring  in  the  filtered 

result .  Presently  it  appears  that  the  most  favorable 

adaptive  approach  is  a  multiple  model  doubly  stoc- 

-JTROD UCTI0 tj  hastic  procedure  (Ref.  2)/which  combines  a  set  of 

local  Markov  models  of  the  field  itself  with  a 
lower  level  2-D  Markov  chain  for  the  model  parameter 
transitions.  This  development  differs  from  pre¬ 
viously  published  procedures  in  that  it  takes  int.o 
account  the  need  for  transition  probabilities  be¬ 
tween  local  models,  rapid  edge  detection  to  prevent 
blurring,  and  near  optimal  local  recursive  estima¬ 
tion  for  general  AR  models. 


In  the  past  several  years  various  methods  have 
been  proposed  for  the  recursive  estimation  of  images 
and  other  two-dimensional  (2-D)  data.  .Most  of  these 
methods  have  been  tied  to  first  order  or  separable 
models  that  do  not  match  image  statistics  very  well. 
Also  most  of  these  methods  are  restricted  to  the 
homogeneous  case  where  the  model  coefficients  are 
constant . 


2-D  RECURSIVE  ESTIMATION 


Since  many  2-D  data  fields ,  including  images , 
are  markedly  non-homogenecus ,  there  is  a  need  for 
general  recursive  estimation  procedures  which  can 
take  this  property  into  account.  The  2-D  Kalman 
filtering  methods  can  theoretically  take  this  into 
account,  however  the  practical  problem  remains  of 
obtaining  the  spatially  varying  model  coefficients 
Thus  recent  work  has  been  concerned  with  the 
development  of  various  adaptive  estimators  which 
will  provide  a  practical  means  of  estimating  the 
model  coefficients.  A  common  property  of  these 
aiaptive  estimators  is  that  they  separate  ir.tc  two 
parts:  one  part  estimates  the  model  coefficients 
based  directly  on  the  noisy  iata  and  a  second  part 
filters  the  data  using  the  estimated  model  coeffi¬ 
cients. 


Early  attempts  to  achieve  a  truely  recursive  g 

2-D  Kalman  filter  were  of  only  limited  success  due  3 

to  both  the  difficulty  in  establishing  a  suitable  & 

2-D  model  and  also  the  high  dimension  of  the  re-  yj 

suiting  state  vector.  In  (Ref.  15)  Voods  and 
Radevan  presented  two  new  algorithms  which,  to  a  ^ 
large  extent,  overcame  the  computational  problems  m 
which  had  precluded  the  use  of  2-D  Kalman  like  pro-  to 
cessors.  Both  vector  and  scalar  scanning  methcds  ^  j 
were  considered,  but  emphasis  was  placed  on  the  0*  < 

scalar  scan  because  it  leads  to  processors  which  OT  : 
are  recursive  in  both  dimensions,  i.e.  2-D  recursive  a;  S 
filters.  These  Kalman  based  algorithms  allow  the  **  6 
use  of  space-variant  models  which  can  provide  a  bet¬ 
ter  match  to  local  source  statistics  leading  to  a 
greater  noise  reduction  with  less  signal  distortion. 


Our  original  approach  was  based  on  the  use  of 
continuously  updated  mcdel  parameter  estimates  to 


Other  approaches  to  recursive  image  estimation 
include  the  work  of  Strintzis  (Ref.  U),  Jain  (Ref. 

5)  and  Murphy  and,  Silverman  (Ref.  6).  In  (Ref.  U), 
Strintzis  presents  an, APMA  modeling  approach  to  re¬ 
cursive  processing  coupled  with  an  interesting 
additive  rather  than  multiplicative  based  spectral 
decomposition.  The  addition  of  a  moving  average 
to  the  AR  model  of  Kalman  filtering  can  provide 

a  tf*  U  ^  SS&2T*! 
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some  modeling  advantages,  however  in  (Ref.  L)  the 
developed  model  is  for  the  correlation  functions , 
not  the  field  itself,  thus  simulation  is  not  pos¬ 
sible.  Further  only  the  steady  3tate  is  considered, 
where  the  image  filtering  problem  can  be  embedded 
in  the  class  of  cyclic  processes.  Unfortunately 
this  embedding  cannot  be  done  if  non-homogeneous 
image  models  are  considered  or  if  random  boundary 
conditions  are  explicitly  treated,  as  is  done  in 
(Ref.  7). 

In  (Ref.  5),  Jain  considers  both  implicit  and 
explicit  partial  difference  equation  models  for 
images.  The  models  are  compared  and  evaluated 
for  mean-square  improvement  on  representative 
images.  This  method  is  not  suitable  for  space- 
variant  or  non-homogeneous  image  models.  In  (Ref. 

6)  Murphy  and  Silverman  consider  a  vector  scanning 
approach  to  image  restoration.  A  constrained  or 
reduced  update  vector  processor  i3  developed  which 
nicely  compliments  the  scalar  reduced  update  filter 
of  Woods  and  Radevan  (Ref.  3).-  The  problem  of 
image  restoration  rather  than  simple  noise  filter¬ 
ing  is  considered  in  (Ref.  6).  In  (Raf.  7),  , 

Woods  and  Ingle  extend  the  reduced  update  filter 
to  the  deconvolution  case. 

ADAPTIVE  ESTIMATION 

Development  of  a  procedure  for  online  parameter 
identification  requires  either  a  recursive  type 
algorithm  in  which  each  estimate  is  a  simply  evalu¬ 
ated  function  of  the  previous  estimate  and  the 
current  measurement,  or  a  non-recursive  algorithm 
which  can  be  executed  in  a  time  frame  significantly 
less  than  a  single  pixel  scanning  period. 


Below  we  summarise  the  development  of  a  corres¬ 
ponding  multiple  model  procedure  for  the  estimation 
of  images  and  other  2-D  fields.  For  more  details 
the  reader  is  referred  to  (Refs.  2,  17,  1 8). 

MULTIPLE  MODEL  RECURSIVE  ESTIMATION 

Following  the  notation  in  (Ref.  3),  the  dynamic 
model  of  the  image  formation  is  given  as, 

s(m,n)  =  l  cij  s(m-i ,  n-J )  ♦  v(m,a)  (l) 

R«+ 

where  the  coefficients  are  the  model  parameters, 

R-.  is  the  non-symmetric  halfplane  support  of  the 

model,  and  w(m,n)  is  a  white  Gaussian  noise  having 
zero  mean  and  variance  a£.  The  scalar  observation 

model  is  given  by: 

r(m,n)  =  s(m,.i>*v(m,n)  (2) 

where,  v(m,n)  is  a  white  Gaussian  noise  with  zero 
mean  and  known  variance  Then  using  the  global 

state  vector  s_(m,n)  as  shown  in  Fig.  1,  the  vector 
equations  for  Eq.  1  and  Eq.  2  can  be  written  as: 

s_(m,n)  =  C  s_(m-l,n)  +  w(m,n)  (3) 

T 

r(a,n)  *  h  s(m,n)  +  v(m,n)  (l») 

where,  C  is  the  system  propagation  matrix  determined 
by  {c  UT  and  the  ordering  of  the  global  state  vector, 

h  «  [1,  0,  0,  .  .  . ,  0]T  (5) 


With  regard  to  developing  recursive  algorithms 
suitable  for  use  on  2-D  fields,  some  type  of  approxi¬ 
mation  is  desirable  because  of  the  nonlinear  nature 
of  simultaneously  estimating  both  states  and  para¬ 
meters.  For  example,  a  recursive  least  squares  or 
weighted  least  squares  algorithm  (Refs.  8,  9)  might 
be  applied  directly  to  the  measurements  in  order 
to  obtain  parameter  estimates  which  can  then  be 
used  for  Kalman  gain  adjustment.  However,  because 
these  procedures  give  biased  parameter  estimates , 

Sen  and  Sinha  (Ref.  10)  have  proposed  a  revised 
least  squares  procedure  which  results  in  approxi¬ 
mately  uncorrelated  residuals.  Alternate  least 
squares  type  methods  for  reducing  the  bias  include 
algorithms  based  on  correlating  the  measurements 
(Ref.  11)  and  those  based  on  correlating  the 
estimates  (Ref.  12). 


and 

iT» 

w(m,n)  *  (v(m,n),  0,  0,  .  .  .,  G]'  .  (6) 

Let  s^m.n)  be  the  local  state  vector  as  shown 
in  Fig.  1,  and  let  r^m.n)  and  v^tm.n)  be  the  ob¬ 
servation  vector  and  white  Gaussian  noise  vector 
respectively  having  the  same  support  as  that  of  the 
local  state  vector  s^(m,n).  Note  that  the  scalars 

r(m,n)  and  v(m,n)  in  eq.  2  are  the  leading  elements 
of  the  vectors  r, (m,n)  and  v^(ra,n)  respectively. 

Then  from  eq.  2, 

r,  (m,n)  *  ^(n^n)  +  v,  (m,n)  (7) 

BASIC  ADAPTIVE  ALGORITHM 


An  alternative  to  the  above  estimation  based 
procedures  which  appears  to  be  very  powerful  in 
view  of  recent  results  in  adaptive  control  (Ref3. 
13,  l1*,  15),  i3  based  on  the  use  of  a  bank  of 
multiple  models  as  per  Lainiotis  Partitioning 
Theorem  (Ref.  16).  Basically  the  Theorem  says  that 
under  certain  assumptions  the  nonlinear  adaptive 
filter  can  be  decomposed  Into  two  parts:  a  linear 
nonadaptive  part  consisting  of  a  bank  of  Kalman 
filters,  each  conditioned  on  a  predefined  model, 
and  a  nonlinear  adaptive  part  which  evaluates  the 
a  posteriori  model  probabilities  and  forms  the  ap¬ 
propriate  state  estimate  as  a  weighted  sum  of  the 
filter  outputs. 


Now  we  make  the  following  hypotheses  to  arrive 
at  a  basic  multiple  model  description  of  an  image: 


1.  At  each  pixel  (m,n),  there  are  L  a-priori  known 
classes  {9j}^^  of  local  state  vectors  distinguished 
by  their  direction  of  predominant  correlation. 

2.  The  probability  distribution  of  the  classes 
(p(9j)}jsl  is  known  a-priori. 

3.  The  conditional  distribution  of  the  local  state 
vector  of  a  given  class  is  Gaussian  i.e.. 


3 


>/2: 


P  (SjJ&j)  %  (2ir)#Vi  |Rj  |  *1/2exp(-l/2  s^Rj  s^) 

(8) 


where  *  dimension  of  s^ . 


where  L(m,n)  is  the  2-D  Markov  chain  field  with 
local  state  support  consisting  of  a  finite  extent 
NSKP  region  similar  to  the  local  state  support  of 
Eq.  1.  Thus  given  {L(m,n)},  the  resulting  signal 
model  is  space-variant  Gaussian. 


Hence  in  light  of  these  hypotheses,  the  dynamical 
equation  becomes 

s_(m,n)  =•  £(9)  s_(m-l,n)  +  w(m,n)  (9) 

where  9  takes  on  one  of  L  values  at  each  (m,n). 

How  a  bank  of  reduced  update  Kalman  filters 
running  in  parallel  and  each  designed  based  upon 
the  statistics  of  one  of  the  L  models  can  be  de¬ 
signed  (Ref.  3).  Using  the  results  developed  in 
{Ref.  17),  the  decision  logic  is  :  Select  model 
8  iff: 

1/2  rf(R  *  ol  IT1?.  +c,  <1/2  r?(R  +  dfll^r.+c.  , 
j  v  — I  J  —  *“ 1  —  i  v  — 1  1  * 

i*J  (10) 

where  c,  •  l/2in|R  +a^I|-inP(9)  (11) 

are  constants  which  can  be  precomputed. 

The  above  adaptive  approach  to  the  use  of 
multiple  models  may  be  improved  upon  by  introducing 
a  spatiaJ.  doubly  stochastic  Gaussian  model  as  dis¬ 
cussed  next. 

DOUBLY  STOCHASTIC  GAUSS I AH  5STIMATI0H 


To  introduce  a  doubly  stochastic  model,  it  is 
necessary  to  have  a  description  for  the  underlying 
process  which  determines  the  elementary  model  to  be 
used  at  each  pixel.  Such  a  description  is  the  2-D 
Markov  chain  which  generates  a  discrete  valued 
random  field  L{m,n)e[3,]  specifying  which  of  the  L 

models  is  to  be  used  at  pixel  (m,n). 


The  simplest  example  of  sufficient  generality 
would  be  a  2-D  Markov  chain  with  local  state  as 
shown  in  Fig.  2.  Such  a  local  state  would  require 
an  L4xL  stochastic  matrix  to  specify  its  transition 
probabilities.  For  the  experimental  cases  to  fol¬ 
low,  L=5  was  found  to  work  rather  well.  In  this 
o a3e  the  transition  matrix  is  625x5.  Each  row  of 
the  transition  matrix  indicates  the  conditional 
probability  of  going  from  each  of  the  IT  possible 
present  local  states  at  !m-l,n)  to  the  conditionally 
possible  local  states  at  (m,n).  The  2-D  Markov 
chain  is  thus  specified  by  giving  the  support  of  the 
local  state  and  its  transition  matrix  P. 


Using  the  2-D  Markov  chain  to  specify  an  'under¬ 
lying  structure,  one  may  generate  a  doubly  stochastic 

l 

model  by  using  L  different  sets  of  parameters  (e,^} 

to  generate  the  random  signal  field.  The  resulting 
signal  model  would  become 

l(m,n)  .  ,  ,»  .  , 

Cj  j  3(m-i ,n-j )  +  u(m,n) 

9+ 


The  filtering  for  such  compound  models  can  use¬ 
fully  be  thought  of  as  a  two  step  process;  at  a 
given  pixel  one  first  estimates  the  local  state  of 
the  underlying  Markov  chain,  then  one  chooses  the 
most  likely  model  to  do  the  estimation  of  the  higher 
level  random  field.  This  will  lead  to  a  processor 
similar  to  the  basic  adaptive  algorithm  but  with 
the  added  advantage  that  the  local  state  of  the 
underlying  Markov  chain  can  alter  the  likelihood  of 
the  various  model  transitions  and  hopefully  improve 
upon  model  switching  decisions. 

Accordingly  assume  that  s  and  5.  are  known  for 
the  past  at  pixel  (m,n).  The  probability  of  L(m,n), 
conditioned  on  the  noisy  observation  of  the  local 
state  r^Cm.ni.can  be  written  as 

p[2(m,n)  jr^m.n) ,  past  .1]  * 
pCr^m.n)  |L(m,n) ,  past  2,]P[L(m,n)  |  i,  (m-1  ,n)  ] 

l  p[r  (m,n)|£,  past  2}PU|L  (m-l,n)] 

1= l  1  -1 


where  use  has  been  made  of  the  conditional  probabi¬ 
lity  of  the  2-D  Markov  chain.  The  function  P( • | • ) 
then  is  Just  given  by  the  corresponding  row  of  the 
probability  transition  matrix  of  the  chain.  The 
conditional  probability  density  of  r^(m,n)  is 

Gaussian  with  mean  cero  and  covariance  determined 
by  the  past  and  present  chain  data  i.  If  we  assume 
all  nearby  points  are  in  the  same  state  1,  then  we 
may  approximate  the  conditional  density  of  r.  as 
2 

H{0,R^+Cv  I_) .  If  we  replace  3.^(m-l,n)  by  its 
estimate  i^Cm-l.n),  we  obtain  the  Jollowing 
decision  directed  rule: 


Choose  S.  such  that 


i/s  gthyl 


Ij'^r^+CjJltm-l.n))  (lit) 


where 

c^)  6  1/2  tn|Rj-K?  II  -  in  P[2.  |  ^ (m-1  ,n )  I  (15) 


This  decision  rule  can  be  compared  with  Eq. 

10  and  11,  the  decision  rule  for  basic  algorithm. 

The  difference  is  the  replacement  of  the  a-priori 
probabilities  P(9^)  by  the  conditional  probabili¬ 
ties  of  the  Majkcv  chain  given  the  decision-directed 
present  state  £^(m-l,n).  If  the  recent  decisions 

have  been  good,  this  has  the  effect  of  improving 
the  probability  of  the  correct  decision  on  the  £ 
value  at  the  present  pixel. 


s(m,n) 


(12) 


m 
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RESULTS 

The  results  of  statistical  Measurements  on 
sections  of  real  images  show  that  a  good  description 
of  the  images  can  be  obtained  with  as  few  as  five 
models.  Four  of  the  models  correspond  to  the  pre¬ 
dominant  correlation  links  in  four  directions  making 
angles  of  0°,  45° >  90°  and  135°  to  the  horizontal, 
for  the  images  sampled  and  processed  as  rectangular 
arrays.  The  fifth  model  describes  the  parts  of 
the  image  with  isotropic  structure.  To  test  our 
algorithm  we  used  a  geometric  image  having  the 
above  mentioned  four  predominant  edges  and  an 
isotropic  structure. 

This  geometric- image  was  a  (128x129)  density 
domain  image,  quantised  to  8  bits  (256  gray-levels), 
(Fig.  3a).  For  simulation  purposes  a  signal-to- 
noi3e  ratio  of  3dB  was  set  up  with  signal  variance 
of  'unity  and  mean  of  zero  (Fig.  3b).  Second  order 
steady  3tate  filters  were  designed  for  all  models 
by  using  the  reduced  update  error  covariance  equa¬ 
tions.  Fig.  3c^ shows  the  estimate  obtained  by  the 
doubly  stochastic  Gaussian  algorithm.  The  signal- 
to-noise  ratio  of  Fig.  3c  is  13.96  dB,  thus  the 
improvement  with  respect  to  the  3dB  input  image  is 
1=10.86  dB. 

We  also  processed  a  face  image  (Fig.  4a)  at 
3ignal-to-noise  ratio  3  dB.  The  noisy  density 
domain  image  i3  shown  in  Fig.  4b.  The  five  models 
were  determined  from  the  noise-free  lata  of  Fig.  4a 
using  a  mask  determined  by  using  the  decision  logic 
of  Sq.  10  based  on  the  geometrical  image  models. 

Thus  this  is  the  result  of  the  first  step  in  a  pos¬ 
sible  iterative  procedure  which  can  tune  the  edge 
and  isotropic  models  to  a  particular  image.  The 
result  is  3een  in  Fig.  4c  which  has  a  signal-to- 
noise  rat'o  of  12.4  dB  which  is  equivalent  to  an 
improvement  1=9.4  dB. 

Example  comparisons  with  constant  coefficient, 
linear  estimate  results  are  contained  in  references 
(2,  IT,  18)  and  the  reader  is  referred  there  to  see 
those  pictures. 

CONCLUSIONS 

The  recently  developed  doubly-stochastic 
Gaussian  estimator  has  been  summarised  and  presented 
as  a  logical  outgrowth  of  a  basic  adaptive  algorithm 
based  on  the  Partitioning  Theorem  and  the  2-D  re¬ 
duced  update  filters. 

Experimental  results  were  presented  U3ing  the 
new  estimator  for  two  test  images.  The  subjective 
and  numerical  results  demonstrate  the  utility  of 
the  adaptive  approach. 
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